function G = gaussCoeGmat(coe,mu)
% RTN frame

a = coe(1);
e = coe(2);
incl = coe(3);
W = coe(4);
w = coe(5);
f = coe(6);
% theta = f + w;

p = a*(1-e^2);
r = p/(1+e*cos(f));
h = sqrt(mu*p);
eta = sqrt(1-e^2);
% n = sqrt(mu/a^3);
% v = sqrt(2*mu/r-mu/a);

G = zeros(6, 3);

% if RTN == 1
ga_r = 1*2*a^2/h*e*sin(f);
ga_u = 1*2*a^2/h*(1+e*cos(f));

ge_r = sqrt(p/mu)*sin(f);
ge_u = sqrt(p/mu)*((1+r/p)*cos(f) + e*r/p);

gi_h = r*cos(f + w)/h;
gW_h =  r*sin(f + w)/h/sin(incl);

gw_r = -1/h/e * p*cos(f);
gw_u =  1/h/e *(p+r)*sin(f);
gw_h = -r/h*sin(f + w)/tan(incl);

%     gTA_r =  1/h/e*p*cos(f);
%     gTA_u =  -1/h/e*(p + r)*sin(f);
gMA_r =  eta/h/e*(p*cos(f)-2*r*e);
gMA_u = -eta/h/e*(p+r)*sin(f);

G(1, 1) = ga_r;
G(1, 2) = ga_u;

G(2, 1) = ge_r;
G(2, 2) = ge_u;

G(3, 3) = gi_h;

G(4, 3) = gW_h;

G(5, 1) = gw_r;
G(5, 2) = gw_u;
G(5, 3) = gw_h;

G(6, 1) = gMA_r;
G(6, 2) = gMA_u;

% else
%     G = [0                     1*2*a^2*v/mu              0;
%        1/v*(r/a*sin(f))          1/v*(2*(e+cos(f)))       0;
%        0                         0                  r*cos(theta)/(h);
%        0                         0                  r*sin(theta)/(h*sin(incl));
%       -1/e/v*(2*e+r/a*cos(f))    1/e/v*2*sin(f)    -r*sin(theta)*cot(incl)/h;
%        1/e/v*(2*e+r/a*cos(f))    -1/e/v*2*sin(f)     0];
% end

% G = G*Tunit;
